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The effect of signal integration through cis-regulatory modules (CRMs) on synchronization and 
clustering of populations of two-component genetic oscillators coupled by quorum sensing is in 
detail investigated. We find that the CRMs play an important role in achieving synchronization and 
clustering. For this, we investigate 6 possible cis-regulatory input functions (CRIFs) with AND, OR, 
ANDN, ORN, XOR, and EQU types of responses in two possible kinds of cell-to-cell communications: 
activator-regulated communication (i.e., the autoinducer regulates the activator) and repressor- 
regulated communication (i.e., the autoinducer regulates the repressor). Both theoretical analysis 
and numerical simulation show that different CRMs drive fundamentally different cellular patterns, 
such as complete synchronization, various cluster-balanced states and several cluster-nonbalanced 
states. 

PACS numbers: 87.18.-h, 05.45.Xt, 87.16.Yc 



I. INTRODUCTION 

Decoupling simple networks from their native yet often complex biological settings can lead 
to valuable information regarding evolutionary design principles. This motivates the design and 
construction of synthetic genetic networks resembling submodules of natural circuitry in vivo, 
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which in turn lead to the construction of devices and softwares capable of performing elaborate 
functions in living cells [1]. Due to recent advances in bioengineering technology, several prototype 
synthetic genetic motifs, such as logic gates y, 3|, toggle switches 4, 5|, and oscillators 6, 7, 8| 
have been successfully constructed. These simple architectures are thought of as essential modules 
in living organisms, and based on them, complementary approaches have been developed to explore 



the relationship between the structure and function of more complex genetic circuits [9|, [1 

A natural step in the design of artificial gene networks would be to include a mechanism of 
intercell coupling that would globally enhance, given that cells are frequently subject to chem- 
ical signals from neighboring cells, the oscillating response of the system. The most common 
communication mechanism with such a function is quorum sensing, the ability of bacteria to 
communicate with each other through signaling molecules that are released into the cellular envi- 
ronment. Quorum sensing has lead to progra mm ed p opulation control in a bacterial population 
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29| . Through such a mechanism. 



the ability of cells to communicate to one another allows them to coordinate the behavior of the en- 
tire community, where gene expression is regulated in response to the local cell population density 
30|. A well-defined example of coordinated global behavior in bacteria is a population of genetic 
relaxation oscillators coupled to a quorum sensing apparatus, which can achieve synchronization 
through the so-called "fast threshold modulation" mechanism 



18[ . Coupling, however, can be 
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devised in different ways, e.g., attractive or repulsive cell-to-cell communication 
synthetic systems. Different couplings would lead to different dynamic patterns, such as synchro- 
nization, clustering and multistability 



28| 



Why is there such a difference in cellular patterns when different types of cellular communication 
are employed? Actually, biological functions appearing as collective behaviors may arise from a 
particular module that integrates intracellular and extracellular signals. Such a module is now 
known as cis-regulatory module (CRM), which contains a cluster of binding sites for transcription 
factors (TFs) and determine the place and timing of gene action within the network, e.g., the 
CRM in the sea urchin embryo can control not only static spatial assignment in development 



but also dynamic regulatory patterning [31'|. TFs are often integrated in a combinatorial logic 



manner, and moreover such a combination may take different schemes 
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leading to different CRMs. In fact, from views of evolutionism, CRMs are changeable, e.g., cis- 



regulatory mutations 38| , and such a mutation constitutes an important part of the genetic basis 
for adaptation. However, how different CRMs affect collective behaviors across ensembles of genetic 
oscillators with cell-to-cell communication remains to be fully explored. 

In this paper, we investigate this question in detail and find that CRMs play a significant role 
in the mode of dynamic patterns at the cellular population level, e.g., the CRMs can drive fun- 
damentally different cellular patterns such as synchronization and clustering. We first design and 
construct a multicellular network with a CRM, using a variant of the synthetic genetic relaxation 
oscillator developed in E. coli ,8] and utilizing quorum sensing to communicate between cells. Since 
different CRMs due to cis-regulatory mutations 38| lead to different types of cis-regulatory input 
functions (CRIFs) such as AND, OR, ANDN, ORN, XOR, EQU, we then investigate the effects of 
these different CRIFs on cellular patterns to support our conclusion. We emphasize that since the 
proposed genetic relaxation oscillator is composed of interacting positive and negative feedback 
loops, and this circuit to polo gy is common in genetic oscillators such as cell cycle and circadian 
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45l | , our conclusion on how CRMs influence the dynamics of genetic 



circuits with this shared topology will be of general relevance to a wide range of cellular processes. 



II. MATHEMATICAL MODEL AND THEORETICAL ANALYSIS 



A. Model 



First, we report our design on a network of coupled synthetic genetic relaxation oscillators with 
a CRM, which is schematically shown in Fig. [TJ The core oscillator is a variant of the genetic 
relaxation oscillator proposed in Ref. [8|. In such an oscillator, the activator X (CII) and the 
repressor Y are under the control of different promoters from the A phage virus. In Fig. [TJa), X is 
the autocatalytic portion of the oscillator whereas Y is a protease that degrades X. Both genes x 
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FIG. 1: (Color online) (a) Schematic diagram for the network of genetic relaxation oscillators with a 
cis-regulator module, where the right bidirectional arrow indicates that S can freely diffusive through the 
cellular membrane, (b) Six cis-regulatory constructs for implementations of six different logic functions. 
In (a) and (b), X, Y and Luxl denote the proteins, Pi and P2 represent the promoters. Or stands for 
operator site whereas RNAP for RNA polymerase. We use offset and overlapping boxes to indicate the 
mutual repression and the dashed lines to indicate the cooperative interaction. 



and y are activated by protein X. Such a circuity not only is a useful architecture to understand 
information processing of sirnple oscillators but also appears as a common core motif in biological 



contexts 



3S 



4C 



41 



42 



43 



bacterium Vibrio Fischeri 



44|. In our design, we utilize the quorum-sensing apparatus of the 
3(J| to communicate between cells. This cell-to-cell communication 
system operates by diffusing a small molecule [also called autoinducer (AI)] into the environment. 
Since the communication is implemented by the signal molecule which regulates the activator X, 
we refer to it as activator-regulated communication. When this molecule binds to a regulatory 
protein (LuxR). both it and X bind the regulatory region of gene x or y and combinatorially 



modulate the transcription rate. Many of these combinational effects are performed by a CRM, 
which can function as analogous implementations of logic gates. The corresponding CRM contains 
a cluster of binding sites of two different transcription factors (TFs) that control the activation or 
repression of a gene. These TFs may be either activators enhancing the binding or the activity of 
the RNA polymerase in the cognate promoters, or repressors blocking this binding, or both via the 
mechanism of "regulated recruitment" 



46j . Based on the possible combination of the two TFs, 
the CRM can perform different logic functions with different implementations, as shown in Fig. 
[U^b). Limited by the regulatory structure of the relaxation oscillator (more precisely, the TF X 
serves as activator only), we have six biologically feasible CRM designs: AND, OR, ANDN, ORN, 
XOR and EQU (see Table I). These logic functions have been either described experimentally or 



suggested to occur on the basis of simulations using empirical data [49j . Actually, the prokaryotic 
transcription networks provide a large number of composite logic operators that are implemented 
through more complex natural or simulated regulatory setups. Alternatively, the CRM designs 
can be implemented by introducing mutations at the amino-acid sequences of the TFs and the bp 
sequences of the cis-regulatory regions 3j|. Note that, in our designs, the signaling molecules can 



serve as not only activators but also repressors by the introduction of an alternative promoter 
TABLE L Logic operations for cis-regulatory input functions 
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Then, we define the chemical species in Table IL All biochemical reactions are listed in Fig. 
^a), and some reaction constants are listed in Table III. Assume the fast reactions to be in 
equilibrium, refer to the equilibrium equations shown in Fig. [21Ib), where square brackets stand for 
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FIG. 2: (a) The biochemical reactions are classified as two classes: fast and slow; (b) The equilibrium equa- 
tions for the fast reactions; (c) Six logic operations and their biochemical reactions, where the corresponding 
conservation laws are listed on the bottom. The reaction rates used are experimentally reasonable, refer 

to idl. 



concentrations of species. In fact, the fast reaction equilibrium trick based on quasi-steady state 
approximation approach has been widely applied to reduce the complexity of multiscalc problems 



m 



54| . The conservation laws for DNA binding sites in the regulatory regions are listed on the 



bottom in Fig. [D^c). 





TABLE II: Description of 


species 


in biochemical reactions 


Species Descriptions 


Species Descriptions 
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Define concentrations as our dynamical variables (see Table IV). Using equalities for the fast 
reactions and the conservation laws, we can eliminate fast variables. To that end, we can derive 
expressions of five cis-regulatory input functions (CRIFs) which are listed in Table V, and the rate 
equations which describe the evolution of the concentrations of X, Y, L and S monomers as follows 



dt 



= CRIF ~ 5,yX,Y, - S^X, 



dt ""^ 1 + Xf 



du 1 + mxf 

dSj 

dt 



(1) 



1 + ^/ 
ctsLt. - SgSi + ri{Se - Si), 

where ~ ^ ^iLi (when N cells are considered) in which Q depends on the cell density in a 
nonlinear way. The rescaled parameters are also listed in Table IV. 
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TABLE III: Descriptions and values of raw parameters. 



Dsscriptioiis 


Values 


References 


Diiii6rizcLtion 6C[uilibrium const&nt 


Ki = 1.8 X 10^*M~^, K2 ^ Kj, ^ 2.5 x lO^M"^ 


\A6 501 


R,Ggul3,tory binding const&nt 


7^4 = TV's = _R'6 = A'7 = 5 X lO^M"-^ 
Kii = Ks) ^ 2.5 X 1Q^M~^ 


[461 


Degradation rate of protein 


dx = dy = ln2/10min~^,dL = ln2/0.2min~^ 
ds = ln2/l.lmin~^ 


[511 


Autoinducer synthesis rate 


c = l.lmin^^ 


[52] 


Bulk rate of transcription and translation 


kx ~ ky = kL = kc = 30min^^ 


[46? ] 


Amplified factor of transcription rate 


fx ^fY = fL^ 10, fc = 90, fxc = fcx = 90 


[46] 


The rate of repressor degradation by Y 


Kxy = 2 X lO-^M"^ 


[18] 


Plasmid copy number 


mx ~ 10, my — l,mL — 50 


[46] 


Concentration of LuxR 


[LuxR] = 1 X lO^^M^^ 


[18] 


Other parameters 


nxi^x[D^'^][P] = 8 X 10"*Mmin"\nx = ny = = 1 


18] 



TABLE IV: Rescaled variables and rescaled parameters for models. We assume K2 = K3, K4 = K5 — Kg — 
K7, Ka = K9, Kx = Ky = Kl = Kc, [D^'^] = [D^'^] = [D^T] = [D^T], and nx = ny = ul = ns = e 
for rescaling. 

Rescaled Variables Rescaled Parameters 

X = {K4Kiy^^[X] n^^fx,fiy = fy,^i = fL,iJ,s = fc, ^^xs = fxcKs/Kr + fcxK<)/Ki 

Y = {K4,KiY''^\Y] Ox = mx,ay = myny /nx,m = mLriL/nx ,t* = nxKx{KiK4)^'''^[D^'^][P] 

L ^ {K4Ki)'^*[L] as ^ cK3{K2Kry^^[LvixR]/((KiK4f^*t'-),t^ t't,5^ ^ dx/T',5y ^ dy/r' 

S ^ iC3(/^2/^7)'/'[LuxR][S] 5i ^ dL/T*,5^y ^ Kxy/{{K^K4f'*T*),5s ^ ds/T*,\^ Ks/Kr + K^/Ki 
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TABLE V: Biochemical reactions and cis-regulatory input functions for relaxation oscillators. 
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B. Analysis 



1. Phase reduction approach 



First, we rewrite the final equation of Eq. (1) as the following symmetric form of coupling 

N 



dSj 
~dt 



a 



,Li - 5sSi - 77(1 - Q)Si + X! ^^^^1 ~ 



(2) 



For convenience, the system composed of both the first three equations of Eq. (1) and the equation 



dt 



asLi - 5aSi - ri{l - Q)Si 



(3) 



is called as auxiliary system, which is assumed to generate a sustained oscillation. Then, we 
perform an analytical study of the entire system in the phase model description, which holds in a 
J. 



weak coupling case 



53]. The main steps are as follows. For convenience, we express the system of 



globally coupled oscillators as 

dxi 
~dt 



1 ^ 



(4) 



where = (X„ F,, L„ 5,)^, / = (Fi, F2, F3, ^4)^ with = CRIF - S^yXiY^ - 6^X,, F2 
^y^+x^ - ^yYi, F3 = ai - SiLi and F4 = agLi - 6sSi - 7/(1 - Q)Si, and p{xi, Xj) 
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(0, 0, 0, rjQ{Sj~Si))'^ . Assume that the uncoupled oscillator has period T. By Kuramoto's theorem, 
for a weakly perturbed system we can obtain the corresponding phase model: 

^ = 1 + . Y.p{x,{^.), x,{^,)\ (5) 

where each Xi{4>i) is the point on the limit cycle having phase 0i, the symbol '-'is the dot product 
of two vectors, and 

Z{(l)i) ^ gv&A(i)i(xi). (6) 

Z{(j)i), a phase response function characterizing the phase advance per unit perturbation, is a 
27r-period function with respect to 0^. To study collective properties of the network, such as 
synchronization and clustering, it is convenient to represent each as (f>i — t + -di with the first 
term capturing the fast free-running natural oscillation d4>i / dt = 1 , and the second term capturing 
the slow network-induced build-up of phase derivation from the natural oscillation. Substituting 
the expression of i?,; into Eq. ^ results in 

d-d I ^ 

= -Z,{t + ^0 . Y.p{x,{t -f di), x,{t + A)), (7) 

The classical method of averaging consists in a near-identity change of variables that transforms 
the system into the form 

N 

dt ^ ' N- 



^ + if2H,,iq^,-cb,), (8) 



where Hij {A(f>) represents the interaction function with respect to the phase difference Acf) = 
between two cells, 

1 r 1 '■"'^ 



H,j - </),0 = - / Z,{t) ■ p{x,{t), xj {t + - ^,))dt =— Z,{9)- p{^j ~(t), + 9)d0 (9) 



which can be calculated numerically 56l|. In what follows, we omit subscripts i and j for conve- 
nience. From H{A(j)), we introduce a function: G{A(j)) ~ H{A(p) — H{—A(j)), to determine the 
mode of coupling. If G{A(p) exhibits a positive slope at Aip = 0, i.e., G"(0) > 0, the coupling is 
phase-attractive; If G'(0) < 0, the coupling is phase- repulsive. Such an approach based on the sign 
of G'(0) that depends generally on the intrinsic dynamics of the uncoupled oscillator and on the 
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23| , especially in the case of complex network 



interaction between the oscillators is more effective than that of directly observing the network 
topology in determining the mode of weak coupling 
architectures. 

According to Tables III and Table IV, we can estimate our system parameter values as follows: 
ax = 10, ay = = 50, a^ = 0.4, bx = 0.5, by = 0.5, (5; = 25, bs = 45, b^y = 5, 10, 

[ly ~ 10, [ii — 10, [is — 9, [ixs = 90, A = 1, ?7 = 10, Q = 0.5. For such a set of values, 
numerical simulation verifies that the term X^jLi vQi^j ~ ^i) affects the timing but not the 
amplitude of the auxiliary system for any > 2, so the above analysis is feasible. In addition, we 
emphas ize that for other different experiments on multicellular systems with the quorum sensing 
11, 



la 
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14 



15 



17|, the differences between the rescaled parameter values are not so large 



that they abolish our conclusions. 



2. Determining the stability of balanced clustering 

Balanced clusters mean that N oscillators are divided into M subgroups of the equal cell number 
with each subgroup being synchronized and with the equal phase difference between neighboring 
subgroups. Here, we employ Okuda's approach |57i] to determine the stability of such clusters 
(see the Appendix of this paper for details). In that method, we need to calculate two kinds 
of eigenvalues: one is associated with intra-cluster fluctuations and the other with inter-cluster 
fluctuations, which are denoted by Xp and Xq respectively, where M < p < N— 1 and < q < M~l 
with Ad being the number of clusters presumptively. For convenience, denote by A^^^ and A^^-* the 
N ~ M same eigenvalues Xp and the maximum of the the real parts of (Af — 1) non-zero eigenvalues 
Xq, respectively. By calculation, we find Xp — 1/^^ X^fcio^ r'(27rfc/M),p M, A/+1, . . . , A^— 1, and 
Xq = l/A/X^fclo^ r'(27rfc/Af)(l - exp(-i27rfcg/M)), g = 0, 1, . . . , M - 1, where T{/^(j)) = iJ(-A0). 
Then, the stability of clusterings can be determined by the signs of A*^^^ and A^^-*. Specifically, the 
clustering is stable if both A*^^^ and A'^^-' are negative, and unstable if A'^^-' is positive. In addition, 
if A^^-* is positive and A^^^ is negative, and further if M = N, the A/-cluster (i.e., the splay state) 
are also stable. 
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In the Sees. Ill and IV, we will numerically study cooperative behaviors of coupled 



relaxation oscillators with different CRMs. In contrast to the previous works 
23, 
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28 
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genetic 



21 



22, 



29|, we will show that different CRMs can drive fundamentally different 



dynamic patterns. 



III. CASE OF TWO COUPLED CELLS: PHASE LOCKING 

Synchronization and clustering of genetic oscillators coupled to quorum sensing result from the 
interplay between the intrinsic properties of the individual cells, the type of cellular communication, 
and the network topology. To gain insights into the rules governing dynamic patterns in complex 
networks of cells, here we investigate the case of two coupled genetic oscillators in detail. 

First, based on Eq. ^ the phase model of two coupled oscillators can be characterized by 

where the phase difference is denoted as A0 — 4>2 — The interplay between the two oscillators 
is often described by the evolution of the phase difference A(/), which is determined solely by the 
odd part of the effective coupling function G{A(j>), i.e., i7(A0) — A(/)). That is, the dynamics 
of A(j) is given by 

The zero points of G{A(j)) are the fixed points of Eq. (fTT|) . These fixed points describe the phase- 
locked states of two coupled cells and their stabilities are determined by the sign of slope of the 
curve G{A(j)) at the zero points: A positive slope means that the corresponding fixed point is 
stable, implying that A(f) nearby the fixed point dynamically converges to the fixed point, whereas 
a negative slope means that the fixed point is unstable, implying that Acf) close to the fixed point 
dynamically diverges. The size of the slope determines the convergence or divergence rate at the 
fixed point. The function G(A0) corresponding to five logic operations AND, OR, ANDN, ORN 
or XOR is shown in Fig. [3)Ja)-(e) respectively (here and below we did not investigate the case of 
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FIG. 3: (Color online) The dependence of the function G on phase difference /\cj> in [0, 2-k\ and the 
distribution of its zero points in the case of two coupled oscillators. In (a)-(d), G has three zero points, two 
of which are unstable and the other is one unstable in the cases of AND and OR whereas one is stable and 
the other two are unstable in the cases of ANDN and ORN. In (e), G corresponding to XOR has 4 stable 
zero points and five unstable zero points. In all the cases, filled circles represent stable points and open 
circles do unstable points, and insets display the corresponding H{A<j)) of the auxiliary system. In (f), 
some typical instantaneous distributions of phases are demonstrated in the five cases of logic operations. 

EQU due to the fact that the EQU destroys the dynamics of the core relaxation oscillator in the 
region of biological reasonable parameters, leading to the loss of sustained oscillation), whereas the 
interaction function H{A(f))s is shown in the insets. 
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(1) AND and OR. The function G(A</«) in Figs. [Si;a)-(b) equates to zero at A(/) = with 
positive slope and at — vr with negative slope. Moreover, the zero point A(/i = is the unique 
stable state of Eq. (fTT|) . Therefore, the phase- model analysis predicts that the phase difference 
of any initial values except A0 ~ tt eventually converges to A0 — 0. This result is also verified 
by integrating the original model with various initial values. A typical snapshot is plotted in the 
upper- left of Fig. El^f). Thus, the analysis together with mmrerical simulation shows that AND 
and OR play a role in stabilizing the in-phase synchronization for two coupled cells. In this case, 
the coupling is phase-attractive. 

(2) ANDN and ORN. Equation pT|) has one unstable state ~ and one stable state A(/i = tt, 
both of which correspond to zero points of the function G, as shown in Fig. [21Ic)-(d). The role of 
ANDN and ORN is to stabilize the antiphase state and prevent the in-phase state. More precisely, 
the integration between the intracellular activator and intercellular signaling repressor in our model 
destabilizes the in-phase synchronization. In this case, the coupling is phase-repulsive. A typical 
snapshot is plotted in the upper-right of Fig. El^f). 

(3) XOR. The function G has nine zero points, four of which, denoted by A(^ = Ai, A2 and 
A</) = TT + Ai, TT + A2, respond to stable states of Eq. pT|) and the other five to unstable states, as 
shown in Fig. . The unstable states form the boundaries for the attraction basins of the stable 
states. The role of XOR is to stabilize four out-of-phase states with phase difference A0 = Ai, A2 
and A0 = vr + Ai, tt + A2 respectively, whereas to destabilize the in-phase and anti-phase state. 
Two typical snapshots are shown in the bottom row of Fig. OJf). 

IV. CASE OF A POPULATION OF CELLS: SYNCHRONIZATION AND 

CLUSTERING 

In this section, we investigate the case of N coupled genetic oscillators {N > 2), focusing 
on two dynamical behaviors, i.e., synchronization and clustering, which are ensemble phenomena 
observed commonly in natural and artificial populations of (possibly weakly) interacting oscillators. 
Synchronization is a cooperative in-phase behavior, which has been the subject of numerous studies 



in physics and biology 



55 



58, 



59 



15 



60|, whereas clustering is a fragmentation of the collective 



behavior in locally synchronized but well separated subgroups, which has been also observed in 



numerous contexts with distinct contributions 
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69( 1 . In what 



follows, we investigate balanced clustering and non-balanced clustering separately for clarity (note 
that synchronization is a particular type of clustering, i.e., 1-cluster). 



A. Balanced clustering 

In the analysis part, we present an approach for determining the stability of balanced clustering. 
Here, we display numerical results for balanced clustering. In particular, we show that CRMs of 
the different structure play different roles in the achieving of collective behaviors. 

(1) AND. Figure IlKa) indicates that 1- and 3-cluster states are stable since both Ai and A2 are 
negative. The instantaneous phase distributions on the unit cycle as shown in Fig. verify the 
coexistence of stable complete synchronization and 3-cluster state. 

(2) OR. In this case, the eigenvalues shown in Fig. UJ^c) indicate that only the complete syn- 
chronization (1-cluster) is stable, which is verified by the numerical simulation shown in Fig. Hl^d). 
The analysis together with the numerical simulation shows that the OR plays a role of stabi- 
lizing complete synchronization, i.e., for any initial conditions for these oscillators, the systems 
consequentially evolve into a stable complete synchronization. 

(3) ANDN. The stability analysis of the eigenvalues shown in Fig. [3Ke) reveals that the network 
of coupled oscillators with the ANDN possesses complex cluster-balanced states, e.g., the stable 
3-, 5- and 8-cluster states. These clustering states are numerically implemented as shown in Fig. 

Hf). 

(4) ORN. We give the results on the stability analysis as shown in Fig. |D(g), which indicate 
that the population of oscillators can give rise to more complex cluster-balanced states than those 
displayed in the case of ANDN, e.g., two additional cluster-balanced states, 9- and 1 1-cluster states 
are found. The instantaneous phase distributions of these clustering states on the unit cycle are 
shown in Fig. Hfh). 
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FIG. 4; (Color online) (Left panel) Eigenvalues associated with intra-cluster fluctuations Ai (blue cycle) and 
the maximal real part of non-zero eigenvalues associated with inter-cluster fluctuations A2 (red square) 
as a function of the number of balanced clusters with flve different logic functions; (Right panel) The 
corresponding instantaneous phase distribution of all possible balanced clusters: (a)-(b) 1- and 3-cluster 
states for AND; (c)-(d) 1-cluster state (complete synchronization) for OR; (e)-(f) 3-, 6- and 8-cluster states 
for ANDN; (h)-(i) 3, 6-, 8-, 9- and 11-cluster states for ORN; (j)-(k) 3-cluster state for XOR. Different 
clustering states appear due to different choices of initial conditions but the cell number is flxed as A'' = 792. 



17 




FIG. 5: (Color online) Temporal evolutions of the concentration of X corresponding to ORN in Fig. [l] 
(a) 3-cluster state; (b) 6-cluster state; (c) 8-cluster state; (d) 9-cluster state; (e) 11-cluster state, where 
each cluster state is indicated by different color or an integer. For each obtained cluster state, numerical 
integration begins from an initial condition close to the corresponding clustering, and plot shown begins 
after allowing a transient time of lO'' units. 

(5) XOR. In this case, the system of coupled oscillators possess only a stable cluster-balanced 
state (3-cluster) that can be seen from the sign of two eigenvalues determining the stability (see 
Fig. mji)). A snapshot of the unique balanced clustering is shown in the Fig. [H^j). 

To display cellular patterns more clearly, we also plot all the time courses of the component X 
in the case of ORN, refer to Fig. [5] These cluster states appearing in the cases of different logic 
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operations indicates that different CRMs can drive fundamentally different cellular patterns. 



B. Non-balanced clustering 

Except for balanced clustering as shown in the previous subsection, we also find non-balanced 
clustering. However, finding all non-balanced clusterings is much more difficult than finding all 
balanced clusterings since the former, one needs to search for all stable regions of initial values 
of coupled systems that lead to stable non-balanced clusterings, and this is even impossible only 
with computer simulation when the cell number is large. Here, we mainly want to show that non- 
balanced clusterings are existent in some cases of five logic operations. By numerical simulation, 
we find that different CRMs can drive different types of non-balanced clusterings except for the 
OR case (since the complete synchronization is globally stable). Taking the cases of CRN and 
XOR as examples, we find several typical non-balanced clusterings which are displayed in Fig. [51 
In the case of ORN, we find 5- and 6-cluster states whereas in the case of XOR, we find a 4-cluster 
state and two different types of 5-cluster states. Note that we did not search out all non-balanced 
clusterings, and other types of non-balanced clusterings except for those found are possible, but 
would depend on the number of oscillators and the choice of initial conditions. 



V. EFFECT OF REWIRING NETWORK ON SYNCHRONIZATION AND 

CLUSTERING 

Biological rhythm results from the interplay between the intrinsic properties of the individual 
cells, the properties of the communication, as well as the network topology. Each property may play 
an important role in shaping the emergent synchronous behavior. Except that different CRMs can 
drive different cellular patterns shown in the above two sections, the rewiring architecture of indi- 
vidual cells also may play a significant role in promoting synchronization or antisynchronization of 
coupled cells, e.g., it has been shown that rewired interaction in a repressilator population with cell- 
cell communication can offer diverse dynamics, such as multistability and clustering 
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FIG. 6: (Color online) Non-balanced clusterings found in the cases of ORN and XOR. (left panel) Temporal 
evolutions of the concentration of X corresponding to ORN and XOR and (right panel) cellular patterns: 
(a) 5-cluster state for ORN; (b) 6-cluster state for ORN; (c) 5-cluster state for XOR; (d) 4-cluster state 
for XOR; (e) 5-cluster state for XOR, where each cluster state is indicated by different color or an integer. 
For each shown non-balanced clustering, numerical integration begins from an initial condition close to the 
corresponding clustering, and plot shown begins after allowing a transient time of 10* units. 

that in the investigated-above models, the signal molecule regulates an activator, thus performing 
an activator-regulated communication. Due to biological background of the core genetic oscillator 
and the quorum sensing, however, the signal molecule can also regulates a repressor, leading to 
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so-called repressor-regulated communication in contrast to activator-regulated communication. In 
this section, we investigate the effect of this rewiring architecture of motifs on synchronization and 
clustering. 

In contrast to the scheme of signal integration in the previous two sections (refer to a simplified 
scheme shown in Fig. [Tl^a)), in what follows we rewire the interaction of the signaling molecule and 



its regulated gene inside the cell 
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27[ , as shown in Fig. [TJb) . More precisely, we let the signaling 
molecule AI and the TF X combinatorially regulate the target gene y instead of gene x. Completely 
similarly, we can derive expressions of 6 possible CRIFs (see Table V), and the dynamical equations 
describing the time evolution of the concentrations of X, Y, L and S monomers in the following 
form: 



4 



= a 



dt i + xf 

dY- 

— i = CRIF - 5yY, 
dt " 



dLi 1 + ^iiX^ 



(12) 



^ = UgLi - S^Si + ri{Se - Si), 
dt 

where the CRIFs are similar to those in the case of activator-regulated communication, refer to 
Table V except that parameters ax and fj,xs are replaced by ay and fiys, respectively. In both 
cases, the settings of parameter values are also the same except for iiyg = 90. The numerical 
results are summarized in Fig. [71 where all balanced clusterings are listed in two cases of activator- 
regulated communication and repressor-regulated communication for comparison. From Fig. [71 
we see that different CRMs also can drive fundamentally different cellular patterns in the case of 
repressor-regulated communication, but the wave patterns are different from those in the case of 
activator-regulated communication (Data for comparison are not shown). In addition, we show 
how the odd part of the interaction function H{A(f>), G(A(/)), in the five logic operations, changes 
with phase difference Acf) e [0, 2tt] in Fig. [51 where stable zero points (symboUed by filled circle) 
and unstable zero points (symboUcd by open circle) are shown. Our results suggest that the 
architecture of biological systems might make them particularly evolvable, namely, simple shuffling 
of finely-tuned network architectures may render new functionalities of networks with feedforward 
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FIG. 7: (Color online) Diflterent balanced clusterings of original and rewired genetic architectures, (a) 
activator-regulated communication; (b) repressor-regulated communication. The corresponding clusterings 
for two cases are listed on the bottom, respectively. Note: Only the same balanced clusterings are shown 
for three logic operations in the case of repressor-regulated communication, but different non-balanced 
clusterings are possible (data are not shown here). 



and feedback. 

The rational design of biological networks and pathways promises to reveal ways of rewiring 
cells for new biological functions or of gaining insights into the behavior of natural systems. Much 
of the work to date has focused on the manipulation of transcriptional and post-transcriptional 
elements to create synthetic gene networks with desired functions 



7C 



71 



72? 1. In contrast, 



our present study provides a possible arsenal for designing and constructing a network of genetic 
oscillators with different cellular behavior, indicating that rationally reprogramming integration of 
two input TFs by changing a CRM to activate a targeted gene could be used to induce transition 
among various cellular patterns towards the corresponding desired functions. In spite of this, we 
expect that understanding how different CRMs render different responses for the coupled genetic 
oscillators with quorum sensing would provide a valuable insight into designing new synthetic 
genetic circuits. 
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FIG. 8: (Color online) The dependence of the function G(A<j!>) on phase difference G(A(/)), where its stable 
zero points (filled circle) and unstable zero points (open circle) are shown. 



VI. CONCLUSION AND DISCUSSION 



Using models of synthetic genetic relaxation oscillators coupled by quorum sensing, we have 
shown both analytically and numerically that different CRMs drive fundamentally different cellular 
patterns, such as synchronization and balanced clustering, and non-balanced clustering, by con- 
sidering two types of communications: activator-regulated communication and repressor-regulated 
communication. Specifically, in the case of two coupled oscillators, we have shown that different 
CRMs have marked influences on characteristics of phase- locking processes, e.g., two oscillators 
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can display in-phase, anti-phase and out-of-phase synchronization with a certain constant phase 
difference, depending on the type of CRM. In the case of (> 2) coupled oscillators with activator- 
regulated communication, there are 1- and 3-balanced clusters for AND, only 1-balanced cluster for 
OR, 3-, 6- and 8-balanced clusters for ANDN, 3-, 6-, 8-, 9- and 11-balanced clusters for ORN, and 
only 3-balanced cluster for XOR, whereas in the case of N coupled cells with repressor-regulated 
communication, there are 1-, 2-, 4- and 5-balanced clusters for AND, 1-, 2- and 3-balanced cluster 
for OR, 5- and 8-balanced clusters for ANDN, ORN and XOR. In addition, some non-balanced 
clusters have been also found. These results would provide a strategy for a network of genetic 
oscillators: the selection of cooperative rhythmic manner, e.g. synchronization and clustering, 
is governed by the nature of the integration of the intracellular signal and the secretion of the 
biochemical signals through which the oscillating cells are globally coupled. In particular, ge- 



netic network architecture found in synchronous circadian clocks 73j might be constrained since 
the complete synchronization independent of initial conditions takes place only in the case of OR 
type of response. In addition, our results would imply that multicellular organisms evolve into 
some functional CRMs for particular goals (e.g., cellular patterns) by performing an elaborate 
computation for input TFs. 

We expect that our findings will stimulate further investigations under a more realistic condition 



involving stochasticity 
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74 



75 



and heterogeneity 19| as specified in the following four points: 
(1) In a stochastic environment, wc should consider the stability of the obtained desired dynamic 
pattern. Theoretically, Golomb et al, have shown that the clustering state is stable on the condition 
that noise intensity is below a critical value 
the extent of phase synchronization 



Tsl. 



6l| . On the other hand, the global noise can enhance 
77l |. but also can destroy the clustering state like in slow 
switching 78l|. Therefore, we should carefully design the CRMs structure in the presence of noise 
to preserve the desired dynamic patterns. 

(2) In our model, a population of identical oscillators communicate with a uniform coupling, but 
it would be of great interest to study the influence of the cellular variability and coupling strength 



heterogeneity on the synchronization and clustering. If heterogeneity is sufficiently small compared 
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to the coupling strength, we can treat the system as identical oscillators, and otherwise, the effect 
of heterogeneity should be considered. In fact, it has been shown that heterogeneous coupling 
strength and element variability can make the occurrence of clustering states possible in networks 
of neural oscillators 



79( 1 . Similarly, in our case, heterogeneity would result in synchronization and 
clustering. 

(3) Our results were obtained under the condition that the intercellular communication is rather 



weak. However, it is likely that coupling is stronger than that considered here 
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801. Therefore, 



it would be of interest to analyze dynamical patterns in the case of strong coupling. In this 



case, other modes of complex behaviors such as multistability 



29( 1 ■ oscillation death 
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231 



2J|, inhomogeneous limit cycle 
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8l( may also appear in 



23( 1 . aperiodic oscillation and chaos 
our models. 

(4) We point out that our results are in general robust to changes in parameter values if they are 
not chosen close to the margin of oscillation of the uncoupled oscillator. For a kind of response (e.g., 
the response of AND type), however, modes of clustering possibly depend on parameter values. 
For example, for a set of parameter values given above, two kinds of clustering modes in the case 
of AND have been found and displayed, but for a different set of parameter values, other kinds 
of clustering modes are possible. In addition, in the case that parameter values are chosen close 
to the margin of oscillation of the uncoupled oscillator, the system can display richer dynamical 
behaviors expecting to be further investigated, but Kuramoto's phase reduction approach cannot 
be used. 

In addition, we point out that many theoretical studies have shown that biological oscillators 
intertwined with positive and negative feedback loops should have the following essential require- 



ments 



First, negative feedback is necessary to carry a reaction network back to the 
'starting point' of its oscillation. Second, the negative feedback signal must be sufficiently delayed 
in time so that the chemical reactions do not settle on a stable steady state. Third, the kinetic 
rate laws of the reaction mechanism must be sufficiently 'nonlinear' to destabilize the steady state. 
Fourth, the reactions that produce and consume the interacting chemical species must occur on 
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appropriate timescales that permit the network to generate oscihations. Facing to the complexity 
of gene regulatory networks, these mathematical insights reveal the true nature of gene relaxation 
oscillators. Our core relaxation oscillator can show sustained and robust oscillation under the 
guarantee of the above theoretical results. Especially, our coupled positive and negative feedback 
biological oscillator models rely on a separation of time scales between the two components to 
create relaxation oscillations, i.e., the activator must have fast dynamics than repressor. To that 
end, we can increase the plasmid copy number concentrations as well as degradation rates of ac- 
tivator, where high degradation rate has artificially been implemented by using peptide sequences 
appended to the protein to make it a target for proteases in the cell {g , 



84| . Therefore, it would be 



possible to experimentally demonstrate our circuit design. It would be much more useful to take 
a hybrid approach in which experiments and modeling can be performed in parallel to advance 
one another. In a cyclic fashion, experiments can be used to inform the designs of mathematical 
models, which can in turn be used to make experimentally testable predictions. 

Finally, ongoing structural, biochemical and cell-based studies have begun to reveal several 
common principles by which protein components are used to specifically transmit and process 
information. Our studies demonstrate that these relatively simple principles can be used to rewire 
signaling behaviors in a process that mimics the evolution of new phenotypic responses. We expect 
that our work would motivate the investigations in areas such as development, where epigenetic 
inheritance leads to a persistent phenotypic alteration in response to transient signals, or in ccU-ccU 
communication systems that coordinate the rich complexity of group behaviors. 
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APPENDIX A: OKUDA'S APPROACH 

In this appendix, we define cluster-balanced states and study their stability. Each cluster contains 
the same number of oscillators. Thus, we restrict our attention mainly on symmetric states. 
Assume that the phase model of N oscillators is governed by 

^ = " + ^Er(</,,-0,), (Ai) 

where i = 1,2, - ■ ■ ,N. Although O can be given any value in a suitable moving coordinate, we 
assume f2 = below without explicitly refer to it. First, we define a symmetric M-cluster state as 
the state in which N/M oscillators belong to each of M clusters. Since no randomness is including 
in the system, all the oscillators in a certain cluster should be located at the same phase. Let 
$fe denote the phase of cluster fc (fc = 0, 1, • • • , M — 1). From the phase equation, we obtain the 
equation for $fe as 

M-l 

M 



^ M-l 

= M E - ^i)- (A2) 

=0 

We seek solutions to this equation in the form 

$,=a;Wt+^, (A3) 

which implies that the phases of the M clusters are equally separated and rotate at a constant 
frequency uj^^\ Substituting it into the above phase equation, we find that the solution of the 
above form exists if 

Next, we analyze the stability of the balanced M— cluster state. Let us put S(j)i = (j)i — ^k (where 
i belongs to cluster k) and express the linearized equation for iJ^j as dS^ = J$, where the vector 
notation 5$ = {5(pi,S(f)2, ■ ■ ■ , Sip^) and NxN matrix J have been used. Without loss of generality. 
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we assume that cluster k consists of the oscillators with kN/M < i < {k + l)N/M. Then , we have 



al-poE -f3iE ■■■ -(3m-iE 
-Pm-iE aI-(3oE ••• -Pm-2E 



(A5) 



y -l3iE -(32E ■■■ al-poE J 
where I is the N/M x N/M unit matrix and E\s& matrix of the same dimension whose components 
are all 1, a and Pk are expressed as 



1 / 



k=0 ^ ^ ^ 



(A6) 



and primes indicate the derivative with respect to the argument. Since J is a cyclic matrix, the 
explicit form of the characteristic equation of J can be obtained as 



M-1 



iA7-j|= n 



9=0 



^M-1 



. k=0 



M-1 



N 



M-1 



(A7) 



(A - a)^-^ l[ix-a+^Yl ^ke''-'"^^ = 0, 



9=0 \ fe=0 

where i = In this way, we obtain A'' eigenvalues of J in the form 

M-1 



M ^ \ M 

fe=0 ^ 



(A8) 



N 



M-1 



fc=0 



1 M-1 , 

M 2^ 



fe=0 



M 



1 - e 



-i2iTkq/M 



) (g = 0,l,--- ,M-1). 
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